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We present a new interatomic potential for water captured in a charge-transfer embedded atom 
method (EAM) framework. The potential accounts for explicit, dynamical charge transfer in atoms 
as a function of the local chemical environment. As an initial test of the charge-transfer EAM 
approach for a molecular system, we have constructed a relatively simple version of the potential 
and examined its ability to model the energetics of small water clusters. The excellent agreement 
between our results and current experimental and higher-level quantum computational data signifies 
a successful first step towards developing a unified charge-transfer potential capable of accurately 
describing the polymorphs, dynamics, and complex thermodynamic behavior of water. 
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I. INTRODUCTION 



Developing accurate interatomic potentials is essen- 
tial for modeling the atomistic behavior of materials in 
diverse chemical and physical environments. Broadly 
speaking, potentials can be classified as empirical or semi- 
empirical. Empirical potentials — ranging from the rel- 
atively simple Lennard-Jones and Buckingham poten- 
tials to more elaborate force fields such as CHARMM* 
or RcaxFF— — are typically parameterized to represent a 
set of specific structural or thermodynamic properties 
of a given molecular system or material. Consequently 
they may have limited success in representing other non- 
parameterized properties, or materials whose atomic or 
molecular constituents lie outside the parametrization 
set. Systems where charge transfer effects are important 
present a particular challenge to empirical approaches. 
The simplest models incorporate fixed formal atomic 
charges with distance-dependent charge transfer switch- 
ing functions^ Others define a charge-dependent func- 
tional form^ often a simple quadratic as in the ES+ 
methodj^ and adjust the charges via chemical potential 
equalization. On the other hand, semi-empirical poten- 
tials guided by quantum mechanics (QM) — for example, 
the embedded-atom (EAM) 7 -^ and modified embedded- 
atom methods (MEAM)^ tight-binding (TB) theory^ 
SCC-DFTB/CHARMM^ diatomics-in-moleculesH and 
empirical valence bond (EVB) approachc o 14 ! 15 — depend 
on potential parameters derived from ab initio calcula- 
tions or experimental data. The success of these semi- 
empirical approaches hinges upon the ability of the model 
to assimilate relevant QM and experimental information 
within a functional form that depends on a relatively 
small set of parameters, and can be translated readily 



into computer code for efficient application to large-scale 
simulation systems. In the case of the TB and EVB 
approaches, the parametcrizations generally require the 
specification of a carefully tailored set of basis wavefunc- 
tions, the estimation of corresponding overlap integrals, 
and on-the-fly Hamiltonian diagonalizations. These ex- 
plicit QM steps significantly complicate the construction 
of the potentials, thus limiting the size and chemical 
diversity of the systems to which they can be applied. 
Similarly, SCC-DFTB/CHARMM and other QM/MM 
methods^ require the definition of appropriate auxiliary 
conditions in order to handle boundary-matching, charge 
polarization, and long-range electrostatic interactions be- 
tween the quantum and classical (molecular mechanics) 
regions of the system. Semi-empirical potentials such 
as EAM and MEAM do not involve explicit QM com- 
ponents, but in their present form cannot account for 
non-perturbativc changes in the charge states of atoms. 
Consequently, they are not expected to accurately model 
such important biophysical and materials problems as po- 
lar systems, electron transport, defect-driven charge po- 
larization, fluctuating valence systems, complex oxides, 
and reactive dynamics. 

In this study, we present the first implementation of 
a novel charge-transfer embedded atom (CT-EAM) po- 
tential aimed at addressing the issues outlined above, 
and apply it to the structure and energetics of netu- 
ral water clusters (H20)„, n — 2, . . . , 20. Importantly, 
the new potential incorporates quantum mechanical in- 
formation in the spirit of TB, EVB, and related ap- 
proaches, while preserving the intuitive features, ease of 
parametrization, and extensibility of the EAM. The po- 
tential is based on a multiscalc framework recently de- 
scribed by two of the authors. 17 ' 18 This framework is for- 
mally based in density functional theory (DFT) , 19 i 20 and 
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couples the electronic and atomistic length scales within 
a self-consistent classical potential. A key feature of the 
potential is its dependence on the redistributed atomic 
electron densities — and by extension, charge transfer — 
which vary with the instantaneous configuration of the 
atoms within the molecule or material. The parameter- 
ized charge distributions are derived from ab initio calcu- 
lations. This 'atom-in-moleculc' perspective and associ- 
ated effective charge are at the heart of the CT-EAM 
model framework. A second important aspect is the 
imposition of self-consistency between the atomic elec- 
tron densities appearing in the two physically distinct — 
embedding and electrostatic — components of the model. 
This requirement is intrinsic to the CT-EAM theory, and 
is in contrast to other charge transfer models, including 
some based on the EAM, where different functional forms 
are assigned to nominally identical electron densities. 

For this initial implementation, we focus on water as a 
paradigmatic small molecule system of immense practical 
importance to biomolecular and materials applications. 
Water also represents an extremely challenging test sys- 
tem for any classical potential due to its strong polar fea- 
tures arising from underlying charge transfer and charge 
polarization, and associated many-body effects^ 

In the following section we briefly review the extensive 
literature on classical potentials for water, with empha- 
sis on previous approaches to the treatment of charge 
transfer. We then review the EAM method and describe 
previous attempts to adapt the model to the study of 
charge-transfer systems. This is followed by a summary 
of the key features of the recently proposed CT-EAM 
as implemented in the present work. In Section III, we 
present our potential parameterization, and in Section 
IV, the results of the model for various water clusters. 
The paper concludes with a summary and discussion of 
future work. 



II. BACKGROUND 
A. Water potentials 

Understanding the thermodynamic and structural 
properties of water is crucial to modeling many biological, 
chemical and physical phenomena. Despite its relevance 
and importance, there are still unanswered questions re- 
garding the properties of water polymorphs and their ex- 
act roles in solution chemistry as well as in biological 
processes. Developing an accurate model capable of si- 
multaneously describing the gas phase, liquid, and solid 
state properties of water has presented enormous chal- 
lenges. Ultimately, a complete model should be capable 
of describing diverse phenomena such as ion solvation^ 
electro- photo-^l and thermo— dissociation of water; 
dynamical properties of the liquid ) 26 i 27 and anomalous 
thermodynamics^ 

There have been many previous attempts to develop 
potentials capable of describing the various phases and 



properties of water, with varying degrees of success. 
Comprehensive reviews are available in Refs. 29U30j | . and 
we will not attempt to review the potentials in detail, but 
rather highlight essential features. Most have concen- 
trated on describing liquid water properties such as the 
temperature-density variation, second virial coefficient, 
diffusivity, radial distribution functions and structure 
functions; others have focused on accurately reproducing 
gas-phase spectrosocopic data<21 Some of the best-known 
potentials are essentially empirical in nature (ST2^- 
SPC^ SPC/E^i TIP3P^ TIP4P^ TIP5P2 7 .), while 
others have used ab initio calculations carried out on 
small water clusters (monomer, dimer) for their param- 
eterizations (MCDHO^ SAPT^ NCC^^ MCY^ 
NEMO^ 3 - CC-pol^) or a combination of ab initio and 
experimental data (POL5^ DIM water—). In almost all 
of the potentials, the parameterizations are carried out 
with the implicit assumption that the basic structural 
unit consists of the water monomer/molecule (notable 
exceptions being Halley et al.^ Corrales^ and Voth et 
a^i 15 i 22 ) Typically, the molecule is represented by a col- 
lection of point charges placed at suitable sites so as to 
yield the correct dipole and higher multipole moments 
for liquid water, as well as the structures of small wa- 
ter clusters in some cases. The total energy of a system 
comprised of water molecules is expressed as a sum of 
coulombic and non-coulombic terms. Simpler potentials 
hold the geometry of the water molecule as well as the 
values of the point charges fixe d 33 ' 35 i 36 ' 37 i 49 while more 
realistic potentials allow OH bond flexibility, modeled as 
harmonic and anharmonic oscillators] 50 ! 51 Further, rigid 
molecule models like that of Dang and Chang^ 2 - ASP^ 3 - 
and NEMO^ 3 - include polarization effects by accounting 
for induced dipole moments at every atom site in a self- 
consistent manner, while potentials like TIP4P-FQ2& use 
an approach similar to ES-f^ (see Section II. B) to ac- 
count for polarization. Other potentials that account 
for polarization effects include MCDHO) 3 ^ which uses a 
three-site model in addition to a negative mobile charge 
corresponding to a polarized electron cloud, the diffuse 
charge pair potential model of Guillot and Guissani^ 
Polarflex;21 based on empirical valence bond theory, and 
TTMf^ 5 - which uses smeared charges and dipoles. 



B. The embedded-atom method and charge- 
dependent extensions 

The EAM formulation and extensions such as MEAM 
have been used to successfully model a wide range of 
condensed phase systems, including fee metals^ binary 
alloys ; 57 ' 58 tin^ group IV elements such as Si^ and 
even organic polymers^ - In the basic method, the total 
cohesive energy of a system is expressed as a function 
of a local electron density, with each atom viewed as an 
impurity embedded in a host consisting of the remaining 
atoms. The host electron gas provides both ion-ion in- 
teractions and a volume-dependent energy component. 7 
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The total energy of the system is written as follows: 

N 

-Beam = ^ E i: (1) 

i=l 

where N is the number of atoms and 

^ = #(Pi) + 5 (2) 

fj is an element-dependent embedding function of the 
effective local electron density Pj(R;) at atomic site i, 
and represents the collective many-body effects of the 
remaining atoms in the host material; 4>ij corresponds 
to a pair potential between interacting atoms i and j. 
The inclusion of the many-body term Fi in the energy 
expression makes the EAM significantly different from 
traditional two-body potentials. Moreover, in contrast 
to earlier approaches that utilized bulk volume correc- 
tions, the EAM volume dependence is local to each 
atom^ corresponding to an effective dependence on local 
coordination^ Equivalcntly — in a tight-binding bond 
picture of the EAM — the embedding energy can be un- 
derstood in terms of a local moment approximation to 
the density of states^ 

In the simplest EAM formulation, parameterized 
isolated-atom electron densities are associated with each 
nucleus, and Pi(Ri) is approximated by the sum of the 
tails of all neighboring atom electron densities at site i: 

ni 

^(R,) ~ X>?(R^). (3) 

3=1 

Here n, is the number of nearest neighbors of atom i, and 
p°j corresponds to the isolated atomic electron density of 
neighbor j. Various refinements of ~p t have been devised 
to account for neutral electron density polarizatio n 10 ' 63 i 64 
as well as the inclusion of alloying effects.— 

As noted in the Introduction, EAM-based potentials 
in their original formulation do not account for explicit 
charge dependence or charge transfer. To address this 
limitation in their models of metal-oxide systems, Stre- 
itz and Mintmiro^ proposed an extension of EAM, ES+, 
in which an electrostatic energy term E es was added to 
£eam- E es is defined by the equation 

N N 

£ ra =^£>fe)|-^^, (4) 

i=i *,j=i 

where Ei lon (qi) represents the ionization energy of an iso- 
lated atom i, Qi is its charge, and Vij is the coulomb 
interaction energy. Following the Rappe and Goddard 
QEq model,— Ei° n (qi) is expressed in terms of atomic 
charge atomic electronegativity Xi> an d atomic hard- 
ness Jj° via a second order Taylor series expansion about 
the isolated neutral atom energy Ei lso (0): 

Ei^fa) = Er{0)+ X %+ l -J°ql (5) 



In addition, Vij is expressed in terms of effective electron 
densities Qi and Qj of atoms i, j as 

Vij = JJ giiMiiMMj^gl) dr dr >, (6) 

where R^ and Rj are the position vectors of the atomic 
nuclei, and qi and qj are the atomic charges. The qi 
include screened nuclear and polarized valence electron 
components. The latter are modeled using a shape func- 
tion fi with a fixed (optimized) parametric form. The 
instantaneous charge on each atom, which varies as a 
function of atomic configuration, is obtained via chemi- 
cal potential equalization, and requires the solution of N 
coupled linear equations involving the and a set of 
charge- and interaction-dependent electronegativities \i- 
The original Streitz-Mintmire formulation was used to 
represent atomic interactions in the Al-0 system, and 
correctly predicted elastic and energetic properties in 
the bulk, as well as surface energies and relaxations, 
with reasonable assignments of ionic charges for the Al 
and O atoms^ It was later used successfully in dy- 
namical simulations of the energetics of vacancies in 7- 
aluminaaS and in studies of the oxidation of aluminum 
nanoclustersp* 7 - again with reasonable values for the com- 
puted ionic charges. However, Zhou al. noted that the 
model could not describe the behavior of the a phase of 
AI2O3 under compression, wherein the computed charges 
oscillated between large unphysical values at short in- 
teratomic spacings^ This behavior was attributed to a 
compensating effect on the part of the EAM component 
of the ES+ potential, whose particular parameterization 
effectively constrained the atoms from approaching too 
closely. To address this problem, and also enable the use 
of alternative EAM paramctcrizations within ES+, Zhou 
et al. developed a variant in which a priori empirical 
charge bounds were imposed on the ions in the electro- 
static component. The resulting CTIP-EAM model suc- 
cessfully described cohesive and surface energies, surface 
oxidation, and thin-film growth of various Al/Zr-oxide 
systems. 

A primary limitation of both ES+ and CTIP-EAM 
is that they assume a quadratic Taylor series expansion 
about the nominal ionic charges and are thus valid only 
for reasonably small fluctuations about these valuesi 17 i 18 
This precludes a non-perturbative description of charge 
transfer in reactive systems, and the significant electron 
density rearrangements that are induced by strong inter- 
molecular interactions. It also prevents a proper descrip- 
tion of the dissociation of interacting atomic and molec- 
ular species, since the imposed quadratic dependence on 
charge does not transition smoothly to the correct linear 
dependence at long rangei 69 i 70 i 71 

A second pressing issue is the lack of self-consistency 
in both ES+ and CTIP-EAM, since the ^(r^Ri) ap- 
pearing in the electrostatic component of these poten- 
tials is regarded as formally distinct — and is parameter- 
ized separately from — the EAM electron density pf. This 
must be regarded as problematic in light of the intrinsic 
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long-range, many-body nature of charge polarization and 
charge transfer. At the electronic level, it is well known 
that subtle interactions in the vicinity of quantum me- 
chanical curve crossings* 7 ^ and the concomitant interplay 
between short and long-range electronic correlations, can 
have a profound effect on the details of chemical bonding. 
Indeed, such effects in water have been recently the focus 
of considerable theoretical and experimental interest ^ 
Both the problem of significant charge polarization as 
well as the self-consistency issue in ES+ have been noted 
previously in the context of alumina^ At the atomistic 
level, these intrinsically quantum mechanical effects must 
be properly reflected in the design of the potential if it 
is to accurately describe charge transfer and reactive dy- 
namics. 



C. Charge-transfer embedded atom method 
potential for water 

A potential that addresses both of these issues within 
a density functional-based multiscale formalism has been 
developed recentl y 17 i 18 The formalism unifies all extant 
embedded-atom models within a common theoretical 
framework, and as an immediate consequence, general- 
izes to a fully-interacting, self-consistent charge-transfer 
embedded-atom potential. This potential is used here as 
the starting point for constructing a new charge-transfer 
potential for water. 

For details, we refer the reader to the original pa- 
pers. Here we summarize the central results. First, we 
note that the formalism automatically imposes the re- 
quirement that p equal g in the embedding and electro- 
static components of the potential, and incorporates a 
proper treatment of the long-range dissociation of inter- 
acting subsystems! 71 ' 72 These constraints together effect 
the crucial balance between short- and long-range elec- 
tronic correlations. 

The general CT-EAM form has been shown to be 
derivable from the exact quantum-chemical atom-in- 
molecule (AIM) 7 - 7 - and diatomics-in-molecule (DIM)^ 
Hamiltonians. In this picture, charge-transfer-dependent 
embedding functions correspond to one-atom AIM terms, 
while pair potentials map onto two-atom DIM termsi^ 7 - 
This reformulation suggests a practical approach to 
the parameterization of CT-EAM molecular potentials 
based on resonance state (diabatic charge state) poten- 
tial curvesi 7 - 3 - Here, we adopt the parameterization per- 
spective of Refs. (T8II71II72I . which emphasizes the charge- 
transfer electron densities as fundamental variables. Re- 
gardless of which approach is chosen, however, two key 
elements of the original theory must be modified: the 
form assumed by the background embedding densities, 
and the total cohesive energy expression. 

Consider the background embedding density first. We 
begin by decomposing the total electron density into a 
sum of AIM components, denoted by p|(r;Rj). Here 
r represents an arbitrary point in space for the elec- 



tronic coordinate, and the p*'s are assumed centered on 
the corresponding atomic nuclei. The p* play an anal- 
ogous role in CT-EAM to the isolated electron densi- 
ties pi in the original EAM. We use 'atom-in-molecule' 
as a general term for referring to the p*, whether de- 
rived from molecules, clusters, or solids. In principle, any 
physically-justified AIM decomposition can be used. 77 ' 78 
The sole requirement is that the decomposition satisfy 
p(r;R) = J2i Pii r > R-i), where R = {Rj} is a collective 
variable representing the instantaneous geometry of all 
atoms in the system. 

The CT-EAM forms for the background embedding 
densities and effective atomic charges qi are obtained as 
appropriate weighted averages of the density difference 

A i (r;R) = p(r;R)-^(r;R i ). (7) 

Aj corresponds to the electron density distribution of the 
medium in which the ith atom is embedded. Let xf be a 
weight function yielding the spatially-averaged quantity 
Gf: 

9f = J AifoRJxfdOdr. (8) 

If Xi is constructed so as to project out p*(r;R^) from 
p(r;R), we obtain a uniform average of the density dif- 
ference between p* and p", which is simply the effective 
charge: 

q t = J (tftaRO-pJtaR,-)) dr. (9) 

qi is the localized zeroth order moment (L = 0) of 
Aj(r;R). Note that this relation formally connects the 
AIM densities and effective charges. This is the mech- 
anism through which CT-EAM imposes its requirement 
on the embedding and electrostatic components of the 
potential, that p = g. 

If instead we take xf equal to a (5-function centered on 
atom i, and utilize the density decomposition of p(r), we 
obtain an expression for the embedding density p* : 

p*(Ri) = /^(rjRJ-p^rjRiJltfCr-RiJdr 
« J [p(r;R)-pt(r;Ri)]<5(r-Ri)dr 
« X>i*(Ry)- (10) 

In the second step, we have approximated pi by p*. This 
is reasonable because for purposes of estimating the em- 
bedding density, the difference between the isolated and 
AIM densities for the atom experiencing the embedding 
is comparatively small. Most contemporary EAM calcu- 
lations already implement a similar approximation: pa- 
rameterized functional forms for the p"'s in Eq. ([3]) arc 
included within the overall potential specification, and 
thus effectively serve as AIM p*'s. 
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In light of Eq. ((TOj), the CT-EAM background embed- 
ding density p* (R») at atom i corresponds to the localized 
infinite moment (L = oo) of Aj(r;R). qi and p,*(Ri) are 
thus closely related, each expressible as a distinct local- 
ized moment of Aj(r; R). 

The second modification of the EAM concerns the co- 
hesive energy expression. The CT-EAM generalization 
of Eqs. |P)-© isii^ 



r Mi Mij 

E = H E ^,MFiMpt,M] + 5 E E 



ij,M®ij,M 



M=l 



jjti M=l 



(11) 



In the embedding term, the index M sums over the 
Mi integer charge states that are to be included in the 
model for the zth atom; in the pair interaction term, M 
sums over all My pairs of included charge states. i^Af, 
&ij,M and ~Pi m are charge-transfer generalizations of the 
conventional EAM quantities. The £l% t M and are 
weighting factors for the particular integer charge states 
or combinations of charge states that are instantaneously 
populated for a given system configuration. 

In order to make practical use of the CT-EAM for- 
mulation, it is necessary to choose the number of charge 
states to be included for each atom type, and also the 
parametric functional forms to be used for the embed- 
ding functions i*i m and pah interactions &ij m- These 
choices arc discussed below in Section III. 



III. CT-EAM POTENTIAL FOR WATER 

An obvious concern with the fixed-charge models is 
that they lack the flexibility to describe phenomena 
where the neutral water molecule is not necessarily the 
fundamental structural unit. Even more sophisticated 
approaches such as MEVB and SAPT associate charges 
with fixed molecular and ionic species (water, hydronium 
ion). Additionally, they rely on the specification of ap- 
propriate quantum-mechanical basis states in order to 
compute dynamical charges. These features make such 
models difficult to generalize to the study of larger and 
more complex water-containing systems where charge 
transfer effects are expected to play a significant role. 
Important examples include the dynamics of solvated 
proteins^ water-silica interactions^ energy transduc- 
tion in molecular motor proteins^ and the electronic 
and magnetic properties of exotic materials^ 

The use of the EAM as the starting point of our 
approach means that our perspective is shifted from 
larger molecular building blocks to a more fine-grained 
picture — exact in DIM — of perturbed atoms embedded 
in a many-body medium, and explicit two-body interac- 
tions. The formal basis for the methodology in density 
functional theory implies that CT-EAM potentials are 
in principle capable of describing arbitrary charge states 
and energetics of the atoms in any given local chemical 
environment. 



Using Eq. (fTTj) as our starting point, we will now de- 
velop an environment-dependent potential that is param- 
eterized to reproduce the ground-state energy and geome- 
try of the water monomer and dimer for select geometries 
of these structures. The parameterization incorporates 
charge transfer information derived from ab initio calcu- 
lations on the hydronium and hydroxyl ions, the neutral 
isolated water molecule, and neutral water dimer. 



A. Environment-dependent atomic charges 

In principle, the CT-EAM potential should be formu- 
lated in terms of AIM electron densities and a relatively 
complete set of atomic charge states, as outlined in the 
previous section. In this first application of the theory, 
however, our aim is to explore the capabilities of the 
framework in the simplest possible implementation. We 
therefore adopt the AIM atomic charge as a surrogate for 
the background density at a given atomic site; ultimately, 
it will be necessary to utilize more detailed approxima- 
tions of the AIM spatial distributions, particularly for 
dynamical simulations. As is clear from the discussion 
surrounding Eqs. (|ilj)- (|10p . this approximation is equiv- 
alent to replacing the electron density distribution by 
its localized zeroth order moment. We also assume two 
charge states per atom, as discussed below. 

Given the functional form for the CT-EAM energy 
(Eq. (jlip ). our first task is to develop appropriate param- 
eterizations for the atomic charges qi. The dataset used 
to fit the AIM charges is computed using standard pop- 
ulation analysis techniques in conjunction with ab initio 
calculations. We have used the ab initio software pack- 
age GAMES S§2, at the unrestricted Hartree-Fock (UHF) 
level and with a fairly high-quality basis set — 6-31G** — 
to obtain the Lowdin atomic charges.— The same basis 
set and level of theory were used in all parameterization 
calculations throughout this work, and all charge and 
potential parameters were varied in order to limit model 
estimation errors to less than 0.015 e and 0.02 eV per 
molecule, respectively. Although electron correlation and 
other effects such as zero-point energy corrections are not 
included in these calculations, Maheshwary et al. have 
performed extensive ab initio calculations on water clus- 
ters using HF/6-31G** and concluded that overall trends 
in the variation of energy with cluster size remained unal- 
tered with further improvements in basis set and level of 
theory— Indeed, comparison of their results with recent 
accurate X3LYP hybrid density functional energies, com- 
puted with a much larger aug-cc-pVTZ(-f) basis setj^l 
reveals a nearly identical pattern of variation in stabi- 
lization energy for the most stable geometry of (H 2 0) n 
as a function of cluster size n. As our model systems 
we choose different geometries of: i) the neutral water 
molecule, ii) the hydronium ion H30 + , iii) the OH~ ion, 
and iv) the water dimer, in order to represent diverse 
coordination environments. It is important to bear in 
mind that these structures and geometries are used here 



to parameterize charge rather than energy. We therefore 
expect electron correlation effects to be less important 
than the quality of the basis set. 

The motivation behind using different model systems 
is to ensure that the resulting interatomic potential is 
sufficiently robust to describe different chemical envi- 
ronments that the oxygen and hydrogen atomic species 
might encounter in various water polymorphs. The par- 
ticular choice of the H 3 + and OH~ ions is based on two 
key considerations. First, they provide a coordination 
environment for the hydrogen and oxygen atoms that is 
distinct from neutral interacting H 2 dimers. Second, 
H.30 + and OH~ are the two primary dissociation prod- 
ucts of water in solution, and thus are essential to de- 
scribing chemical reactions involving water. 

A distinguishing feature of our methodology is the 
identification of local clusters within an instantaneous 
configuration of the system, and the indexing of an 
atom's charge based on the kind of cluster to which 
it belongs. Each local cluster is assigned a charge de- 
pending on the number of atoms in the cluster, with 
the charge partitioned among the cluster atoms in a 
geometry-dependent manner. For each oxygen atom, we 
identify the number of hydrogen atoms within a radius 
that is chosen to be 1.5 A. For example, if two hydrogen 
atoms are in close proximity to an oxygen atom, then the 
cluster (oxygen plus two hydrogen atoms) constitutes a 
neutral H2O cluster, while if the number of hydrogens 
surrounding an oxygen atom is three, then the cluster 
is identified as a H30 + with a net cluster charge of +1. 
The total charge on an identified cluster with Ah hydro- 
gen atoms is thus iVn — 2. Once all clusters have been 
identified, the total cluster charge is partitioned among 
the atoms as a function of their relative positions within 
the cluster. We then account for further charge polar- 
ization and charge transfer between neighboring clusters 
by parameterizing the amount of charge transferred be- 
tween two water monomers (constituting a dimer) as a 
function of the hydrogen-bond distance between the two 
monomers. The final total charge on a given atom con- 
sists of both intra-clustcr and inter-cluster contributions; 
this corresponds to its effective AIM charge. The follow- 
ing section provides further details of the charge param- 
eterization procedure. 



B. Model Clusters: H 2 0, H 3 0+, and OH" 

For the three model clusters, we initially obtained the 
equilibrium geometries as given in Table [U Next, we 
varied the geometries of the three systems to obtain 
the Lowdin atomic charges as a function of system ge- 
ometry as shown in Fig. [T] For the water molecule, 
atomic charges for the three vibrational modes (symmet- 
ric stretch, asymmetric stretch, and bending) were ob- 
tained, and the charges were fitted as a function of the 
two OH bond distances and the intramolecular angle. 
The symmetric and asymmetric stretches and contrac- 



(a) 



(b) 




FIG. 1: Geometry of the three model clusters: (a) H2O; (b) 
OH-; (c) H3O+ 



TABLE I: Equilibrium geometry values for the model clusters. 
Angles in degrees; distances in A. 
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tions varied from 70% to 140% of the equilibrium bond 
length (R eq ) at various values (55%-130%) of the equi- 
librium intramolecular angle (0 eq ). In a similar fashion, 
we obtained charges for symmetric deformations (70% to 
140%) of the OH bonds with the three bond angles fixed 
at the equilibrium value for the HaO + ion as well as the 
Lowdin charges on the O and H atoms for the OH~ anion 
for deformations ranging from 70% to 140% of the equi- 
librium OH bond distance. We then fitted the atomic 
charge variations as a function of the relative positions 
of the respective atoms in the cluster. 

We now present the parameterization equations relat- 
ing the variation in atomic charge with respect to cluster 
geometry. As noted previously, the number of atoms in a 
cluster is defined by a central oxygen and the number of 
hydrogen atoms that lie within a specified radial cutoff 
r cut = 1.5 A. 

Consider a cluster with a central oxygen O and Ah > 
1 hydrogen atoms. Let the position vector of the p th 
hydrogen atom with respect to the central oxygen atom 
O be r p . The charge on the p th hydrogen atom is 
expressed as a function of the positions of all atoms in the 
cluster, specifically, 9 p o s , r p and r s , where s corresponds 
to any of the other Ah — 1 hydrogen atoms in the cluster, 
OpOs is the angle between r p and r s , and r s is the distance 
of the s th hydrogen atom from O. We have 



1h = «i + «2 +«& H >21 



(12) 



where 



2i = E[ a ^°p) e " 2rp + P(6pOs)r P e- r » 

+ c(9 p os) sin 2 9 p0s , (13) 



s = l 



/Yh 



q 2 



2 = yi ( r p - r s )d(9 p os) sin 2 0, 



pOsi 



(14) 



s = l 
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TABLE II: Water monomer charge model parameters I. 

Qi (e) « 2 (e) Pi (e) (3 2 (e/A) 71 (e) 72 (e) ti (A) 
-10.032 9.3087 7.7183 3.0558 -1.4124 -6.7189 0.03 



and 

IVh 

1nh>2 = [«2e" 2rp + /3 2 f p e~'' p + 72e~ rp ] sin 2 6 p0s . 

s = l 

(15) 

1n h >2 i s non-zero when TVh > 2. a (3, c, d are functions 
of $ p os, defined in Eqs. (f!?0")) - (|2"3"|) below; a 2 , (3 2 , and 72 
are constants whose values are given in Table [HI For the 
special case where the identified cluster contains only a 
single hydrogen, the charge on the hydrogen is given by 

q p H = a 1 e- 2r -+f3 1 e-^+ 11 , (16) 

where cei, (3\, and 71 are constants specified in Table ITT1 
To prevent energy discontinuities, we utilize a switch- 
ing function S(t) to modulate the calculated charge on 
hydrogen atom p as a function of r p : 

?H = ?H 5 (r p - t cut ), (17) 

where 

5(t) = ^(l-tanh(t/ti)). (18) 

Here is given by Eq. (fT!?)) or Eq. (|16|) depending on 
the oxygen coordination, i cut = 1.41 A, and t\ is given 
in Table HH The role of S is to asymptotically switch 
q^ from q^ to zero at a radius that is less than the clus- 
ter assignment cutoff r cut (note that t cut < r cut ). This 
prevents energy discontinuities when an H crosses from 
outside to inside the r cut boundary. Depending on the 
number of hydrogen atoms TVh in the cluster, the charge 
qo on the oxygen atom in the cluster is 

go = (JVk-2)-X;&- (19) 
p=l 

The functional forms for a, (3, c and d are given by the 
following equations (parameters are listed in Table IITT| : 

a(9) = a x e 6 + a 2 e e ' 2/i + a 3 9, (20) 
13(0) = hOe- 8 + b 2 6e- 26 + b 3 6 2 , (21) 
c(6) = Cl 6e- e + c 2 6e- 2e + c 3 , (22) 

d(6) = d^e - d 2 ) + d 2 d 3 . (23) 

Fig. [5] compares the actual and predicted charges for 
the oxygen atom in the three different clusters for select 
geometries. The fits are very good in each case. 
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FIG. 2: Oxygen charge as a function of OH distance for a 
symmetric variation at 9 eq in H2O, HsO + , and OH - . 




FIG. 3: Water dimer geometry, illustrating hydrogen bonding. 



C. Charge Transfer between Clusters: Water 
Dimer Atomic Charges 

Water polymorphs are characterized by the formation 
of hydrogen bonds between neighboring water molecules. 
Thus, the environment of any atom in bulk is differ- 
ent than when it is part of an isolated water molecule. 
In order to include the effect of a bulk environment on 
the atomic charges, we considered two water molecules 
(dimer) and parameterized the atomic charges for select 
geometries of the dimer (see Fig. [3]). This was done by 
fixing the geometries of the individual water molecules to 
match the equilibrium isolated water geometry and mov- 
ing the two molecules relative to each other along the line 
of hydrogen bonding thb between the two molecules. 

Using the same level of theory as above, the equilib- 
rium dimer geometry was computed, yielding intramolec- 
ular bond distance and bond angles of 0.950 A and 
105.5°, respectively. The computed intermolecular hy- 
drogen bond distance ?*hb was 2.03 A, with an inter- 
molecular bond angle </? (formed between the two oxy- 
gens and the common hydrogen; see Fig. [5]) of 172.3°. 
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TABLE III: Water monomer charge model parameters II. 



ai (e) a 2 (e) o 3 (e) foi (e) b 2 (e/A) 


b 3 (e/A 2 ) 


ci (e/A 2 ) 


c 2 (e/A 3 ) 


c 3 (e/A 3 ) 


dr (e/A) d 2 (e/A) 


0.8158 -3.2198 0.5725 -14.0058 61.8232 


1.3188 


12.2353 


-15.7797 


-3.5928 


0.4430 0.1154 



Note that the equilibrium intramolecular bond distances 
and angles for each molecule are very similar to those of 
a single water molecule, while the intermolccular angle 
corresponds to a nearly linear configuration. 

We next varied thb from 75% to 130% of its equilib- 
rium value, keeping the geometry of the two molecules 
rigid and fixing ip at its equilibrium value. This resulted 
in a finite intermolccular charge transfer between the two 
molecules, such that the donor molecule became nega- 
tively charged relative to the acceptor molecule as a func- 
tion of the common hydrogen position. Based on these 
results, we defined a net intermolecular charge transfer 
dq (donor — > acceptor) between the two clusters, and 
a partitioning — referred to collectively as {dqi} — of this 
charge transfer among the constituent atoms, dq is pa- 
rameterized as: 

dq = a m er b "" r ™, (24) 

where rjjB is the distance between the donor hydrogen 
and acceptor oxygen. Based on our HF calculations, we 
choose the cutoff for charge transfer between clusters to 
be s cut = 2.5 A (beyond this distance the computed dq 
was effectively zero.) As above (cf. Eq. (flT]) ). the switch- 
ing function S modulates dq so as to ensure energy con- 
tinuity; it guarantees that if thb > s cu t, there is no in- 
termolecular charge transfer: 

dq = dq <S(r H B -r gim ). (25) 

The parameters ai m , 6j m , and r q i m in Eqs. I[24p and (|2"B"1) 
are given in Table [TV] dq is partitioned among the atoms 
as follows: 

dq o donor = _ Q5 d ~ ( 26 ) 
dq ^""- =0.75 dq, (27) 

and 

dtp*° n °* = -0.4 dq, (28) 

where O do „ or , O accepto ,., and H donor represent the donor 
oxygen, acceptor oxygen, and donor hydrogen respec- 
tively. For the acceptor molecule, dq Rao '"' ptor is com- 
puted by partitioning [dq — dq° aceB *' t ° r ] equally among 
the constituent hydrogens in that cluster. Similarly, 
[(dq° donor + dq Hdor "> r ) — dq] is distributed equally among 
the hydrogen atoms (other than H-donor) in the donor 
cluster. The total charge qi on the zth atom is given by 

q t = qf + d qil (29) 



TABLE IV: Water dimer charge model parameters. 



aim (e) 


b im (A 1 ) 


Tqim (A) 


1.5812 


1.8222 


2.35 



where qf is the atomic charge due to intramolecular 
charge transfer, calculated from Eqs. ([12"]) - (fig)) (qf = 
q^ for H, and qo f°r O), and dqi is the additional 
atomic charge acquired via intermolccular charge trans- 
fer (Eqs. (2"4 ]) -(|2"g ]l ). As we shall see shortly, all three 
values — qi, qf, and dqi — are needed for computing the 
total energy in CT-EAM. qf and dqi are the background 
density arguments to distinct charge transfer embedding 
functions (see Section III.D), and is used to compute 
the coulomb pair interaction energy. 

The specification of the CT-EAM charge transfer pa- 
ramctcrizations for both intra- and inter-cluster interac- 
tions is now complete. Note that we have assumed that 
each cluster is defined such that a given H atom belongs 
to only one cluster. In particular, the identification of a 
'hydrogen-bonding' H atom implicitly assumes that the 
H atom belongs to one cluster and is hydrogen-bonded 
to the oxygen of the neighboring cluster. 

There are a number of special cases that may arise; 
these are handled as follows. If a hydrogen belongs to 
more than one cluster, we initially treat the clusters sep- 
arately, account for their cluster charges, and add the re- 
spective contributions for the common hydrogen. Then 
we use the inter-cluster charge transfer function to deter- 
mine the charge transfer between the two clusters in each 
direction, and add the results. That is, for the given pair 
of clusters, we consider both scenarios where one cluster 
acts as a donor and the other as an acceptor and vice- 
versa. If more than one donor hydrogen is shared between 
two clusters, we use the same set of charge transfer equa- 
tions to compute two sets of charge transfers between the 
clusters: there is no coupling between them. Finally, the 
charge transfer between clusters is always mediated by 
the hydrogen atoms, irrespective of the relative distances 
between the corresponding oxygens. 

D. Charge-dependent embedding functions 

We have described two intrinsic types of charge trans- 
fer in the water system — inter- and intra-molecular — and 
presented parameterizations for each. These two types of 
charge transfer make distinct contributions to the energy 
through their respective charge-statc-dependent embed- 
ding functions (c/. Eq. (|lip ). We must now consider 
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how to determine appropriate functional forms for these 
embedding functions. 

In the original EAM formulation, the atomic em- 
bedding functions were determined by numerical fits of 
the energy to configurational reference states along a 
symmetric dilatation curve.— Later, as a key aspect of 
MEAM, Baskcs proposed the use of a universal plnp 
functional form, with the density argument normalized 
to a reference state. Baskes rationalized this form by not- 
ing that it gave the correct coordination dependence be- 
tween bond length and energy (bond-order/bond-length 
correlation) for SiJ^ Indeed, MEAM has since proved 
remarkably robust in applications to chemically-diverse 
materials systems.— i 10 i 57 i 59 This suggests that the same 
form may also work well as an ansatz for the charge- 
transfer embedding functions required here. 

An independent rationale for the p In p form comes 
from recent work on ensemble models of charge trans- 
fer for strongly-interacting subsystems 3^ In a resonance- 
state (microscopic) ensemble picture, the equilibrium 
charge transfer within a larger closed system provides a 
measure of the interaction strength between subsystems. 
In the equivalent thermodynamic ensemble, the charge 
transfer parameter maps onto a non-zero electronic tem- 
perature. This temperature is conjugate to the charge- 
density entropy induced by the electronic polarization 
and charge transfer among constituent subsystems. In- 
terpreting the charge transfer in terms of an effective 
electronic temperature suggests using the information- 
theoretic form of the entropy Pi In pi, where i indexes 
the pure states contributing to the ensemble^,) , to model 
the charge-transfer embedding energies. 

In light of the universal nature of the density functional 
electronic theory underlying CT-EAM, we expect the em- 
bedding functional form to be independent of the nature 
(inter- or intra-molecular) of the charge transfer. We 
therefore adopt the p In p form for all charge-transfer em- 
bedding functions — using charges instead of background 
densities as discussed below. Finally, the ensemble for- 
mulation and information-theoretic interpretation both 
suggest that distinct charge transfer contributions should 
enter additively into the overall energy expression; this 
is consistent with the formal result in Eq. (fTTjl . 



E. Embedding function and pair interaction 
par ameterizat ions 

The paramctcrizations we have chosen for the AIM 
charges imply a choice of Mi = 2 for both H and O, and 



Mi. 



3 for the pair interactions (Eqs. 



3). We 



write the net embedding energy contribution Fi of the 
ith atom in terms of qf and dqi as: 



Ft = A t qf \n(e qf ) + Afdq t ln(e % 



(30) 



where e = 1 has dimensions of e~ 2 . We use the square 
of the charge to ensure a positive argument for the log- 
arithm; the additional factors of two arc absorbed into 



the parameterization via the prefactors. The intra- and 
interatomic charge transfer values qf and dqi are used 
in lieu of the nominal background embedding densities 
'Pi Mi M = 1,2. The justification for this comes from the 
common origin of q,; and ~p* in Eq. (jSJ). We have also 
absorbed the weighting factors fi^M and f2jj,M into our 
paramctcrizations (Eqs. ([3171) ~ <|3"5|l ). 

Further insight into Eq. (j3"0)) can be obtained by re- 
garding the first term as corresponding to conventional 
EAM, with the AIM charges within the water monomer 
playing the role of the embedding electron density. This 
term is associated with first-neighbor, intra-molecular 
charge transfer. The second term is then a CT-EAM 
correction for second nearest-neighbor, inter-molecular 
charge transfer. It is interesting to note in this connec- 
tion that a second-nearest-neighbor MEAM has been pro- 
posed recently, aimed at correcting the structural stabil- 
ity and surface energy ordcrings in certain bec metals^ 

The total pair interaction <&y is given by the sum of 
two terms, a classical electrostatic component, Vij (anal- 
ogous to Streitz and Mintmire's Vij, cf. Eq. and a 
non-coulombic component, fcj (cf. Eq. ((2J)): 



(31) 



Since we have chosen to utilize localized zeroth-order mo- 
ment models of the AIM electron densities, the electro- 
static component of $y consists simply of the classical 
coulombic interaction between AIM charges qi and qj , 



Vij = qiqj/Rij. 



(32) 



These charges are constrained to be identical to those ap- 
pearing in the embedding component of the potential, in 
accordance with the CT-EAM self-consistency require- 
ment. The form of the non-coulombic potential is dic- 
tated by energy fits once the charge-dependent compo- 
nents have been determined. These assume a purely re- 
pulsive Born-Mayer-type form for the homonuclear pair 
interactions, and a linear-exponential form for the OH 
interaction. They are similar to the functional forms uti- 
lized for pair interactions in the original EAM^ and are 
given by: 



»00 = aooe 



-4r r O o 



aoHfoH + bone n,r ° H + 

r OH 



and 



(33) 

■S(ron - r cut ), 
(34) 



fern = 2a H He- 2rorHH S(r H H - H cut ). (35) 



In these expressions, tq = 1 has dimensions of A , 
and S(t) is the switching function defined in Eq. (|18[) . 
</>hh and ipoo are purely repulsive. They damp to zero 
beyond their respective cutoffs r cut and H cut (the latter is 
specified in Table fV| . For consistency, r cut is taken to be 
the same value as used above for determining whether 
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r 

FIG. 4: Pair interaction potentials (in eV) as a function of 
internuclear separation r (in A), for H-H, O-O, and O-H. 

an H atom belongs to a particular cluster. This pre- 
vents non-coulombic H-H and O-H interactions between 
atoms in different clusters, for geometries near equilib- 
rium. Note that </>oh has been designed to be very re- 
pulsive at small O-H separations by including an Tq^ 4 
term; this prevents the appearance of unphysical energy 
minima. The pair potentials are plotted in Fig. |4j The 
unusual "coat-hanger" shape of </>oh is a consequence of 
the fact that the pair potentials are parameterized in 
conjunction with the electrostatic term </>y, as part of an 



overall fit (c/. Eqs. (|3Tj) and (|32|) ). The particular shape 
prevents O-H interactions between neighboring clusters. 

The parameterization of the embedding functions and 
non-coulombic interactions was carried out with respect 
to a set of reference energies. We chose the symmet- 
ric mode of the monomer for three different bond an- 
gles (105.9, 100, and 110°) and the equilibrium geometry 
of the dimer as our reference configurations. Energies 
at each geometric configuration were obtained by sub- 
tracting the isolated atom energies from the total energy 
obtained via ab initio UHF 6-31G** calculations using 
GAMESS. The energy of the isolated oxygen computed us- 
ing this basis set was —74.7839 hartrees and that of the 
isolated hydrogen atom equaled —0.5 hartrees. Aq and 
A^ are constants given in Table fVl 

For an oxygen atom O in a cluster with N hydrogens 
(N > 2), 

N N 

A = -Aeo^I Yl sin2 (^ofe) 

i=i k=j+l 

x exp [ - ^r 2 (r 0j - r ok ) 2 sin 2 (6» j0 fe)] , (36) 



where j and k represent the j and k hydrogens in the 
cluster, and Aeo is defined in Table fVl If N = 1, we set 
A = 2A eo . 

If a hydrogen atom p that belongs to a cluster contain- 
ing the oxygen atom s is involved in hydrogen bonding 
with TVoh oxygens of Noh different neighboring clusters, 
then 



A H = A 



EH 



NOH / 

1 + r? ^ [exp ( - 2[1 + cos(^ ps )] 2 ) 

u=i 



1 — tanh 



t-2 



(37) 



where u is the index corresponding to the neighboring 
clusters, ups is the angle between pu and ps, and rj. ?*/j S , 
£2 and Aeh are defined in Table IVl Here we have again 
invoked a switching function in order to avoid energy 
discontinuities. If a given hydrogen atom is not involved 
in hydrogen bonding, then Ah = Aeh- 

Figure [5] depicts the actual (UHF calculations) 
and model-predicted variation in energy of the water 
monomer as a function of OH distance for the symmetric 
mode at the equilibrium angle. Table I VII gives a com- 
parison of the monomer properties as predicted by our 
potential, UHF calculations, and experiment. We used a 
modified 93 BFGS routine with analytic evaluation of gra- 
dients to determine the minimum energy (equilibrium) 
geometry. 



It is evident from the table as well as from Fig. [5] that 
the energetics and the minimum energy structure of the 
monomer are well reproduced. However, the dipole mo- 
ment of the monomer as predicted by our potential is 
significantly lower than experiment. Of course, there 
is no physical reason to expect the Lowdin charges to 
reproduce the dipole moments computed as proper ex- 
pectation values. Indeed, if we choose instead a defini- 
tion of the atomic charge based on a physical observ- 
able (the dipole moment),— we obtain the following ef- 
fective local (static) and nonlocal (dynamic) contribu- 
tions to the atomic charge on oxygen in the monomer: 
Z* oc = MW/r|e q = -0.541 and Z*, = r dZ{ c (r)/dr\ eq = 
—0.229, where r refers in this case to the OH distance, 
and the derivative is evaluated for the symmetric stretch 
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TABLE V: Energy parameters, aoo, etOH, bhh, and 6oh in eV; coh in units eV • A ; Aeo, Aq, Aeh, and A^ in eV/e; H cu t, 
rhs, and ti in A; n is dimensionless. 





«OH 


HHH 


froH 


COH 


Aeo 


^O 


^EH 


At, 


H cu t 




t 2 


V 


25.0 


3.0111 


25.0 


-2.4053 


2.5 x 10~ 6 


-11.429 


0.0 


4.7621 


-0.5 


2.43 


2.1 


0.1 


0.0505 




-7.0 1 ' ' ' 1 ' 1 

0.7 0.8 0.9 1.0 1.1 1.2 1.3 

OH distance (A) 



FIG. 5: UHF and predicted energies for a water monomer for 
the symmetric mode at the equilibrium angle. 



TABLE VI: Monomer equilibrium properties. 





Predicted" 


UHF" 


Expt. 


R eq (A) 


0.9431 


0.9431 


0.957 6 


8eq (deg) 


105.47 


105.99 


104.52 6 


H (D) 


1.23 


2.19 


1.86 c 


E eq (eV) 


-6.51 


-6.52 





a Present work. 
kRef.jll. 
c Ref. Ha 



mode at fixed, computed equilibrium angle. The total 
Born effective charge Z* — —0.770 is given by the sum of 
the local and nonlocal contributions. This value can be 
compared with i?L6wdin = —0.444. Similar results would 
be expected for the dimer, where the Lowdin value for 
the dipole moment is in fortuitously good agreement with 
experiment. 

Table I VIII contains information about the equilibrium 
properties of the water dimer (i.e., its minimum energy 
configuration properties), with the geometry defined in 
Fig. [6l The binding energy U is obtained by subtract- 
ing the two monomer equilibrium energies from the total 
energy. For comparison, we also include the relevant ex- 
perimental and UHF results. Once again, we are able to 
reproduce the dimer properties reasonably well with our 
potential. 

At this stage it is important to recall that the parame- 
ters in our final energy model have been determined so as 
to represent the energetics of select geometries of the wa- 
ter monomer, and to yield the correct minimum energy 



TABLE VII: Equilibrium properties of the water dimer. 
All distances in A; angles in degrees; binding energy U in 
kcal/mol; fj, in D. 





Predicted" 


UHF" 


Expt. 6 


roiHi 


0.941 


0.948 






0.952 


0.942 






1.937 


2.038 




r 2 H 4 y 3 


0.943 


0.944 




ro 1 o 2 


2.886 


2.98 


2.952 


ln 1 1 U 2 


105.44 


105.91 




Zh 3 o 2 h 4 


105.37 


106.31 




^OiH 2 2 


177.3 


179.27 


174.0 




1.8 


-0.5 


0.0 ±6.0 


1> 


60.9 


56.7 


58.0 ±6.0 


u 


-5.860 


-5.505 


-5.40 ±0.7 


II 


2.30 


2.60 


2.64 



"Present work. 
6 Refs.|9li3 




FIG. 6: Structural parameters defining the optimized water 
dimer structure. 



structure of the water dimer. No energetic information 
for the remaining two model clusters used in the charge 
parameterization step — H30 + and OH - — was included 
in this energy parameterization process. Consequently, 
we should not expect the model in its current form to be 
able to accurately predict the energetics of ionic molec- 
ular species. We therefore focus on assessing the predic- 
tions of the potential for the structure and energetics of 
small neutral water clusters. These results are summa- 
rized in the following section. 
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IV. RESULTS 

In developing a model capable of accurately de- 
scribing water polymorphs, a basic but important 
requirement is the ability to predict the correct 
structure and binding energies of neutral vapor- 
phase water clusters^ There have been numerous 
computational^^S^L^i^i^ ^ 

mental studies^I^^Li ^ ^ .^^^ examining 
various water clusters. It has been shown that small neu- 
tral water clusters have 2D cyclic structures, where each 
molecule serves as both an acceptor and a donor, while 
the larger clusters have 3D structures. This crossover is 
seen for the water hexamer and larger clusters, where the 
3D structures are energetically favored. Some of the pop- 
ular water potentials (MCDHO, TIP5P, TIP4P, POL5, 
Dang and Chang (DC)) have been used to study small 
water clusters with varying degrees of success. In the 
following, we compare our results with these potentials, 
as well as experiments and quantum calculations. We 
pay particular attention to the different structures of the 
hexamer. 



Trimer ring 



\ /20Q 

(a) 6 



Tetramer 



r 



(c) 



Trimer chain 



(b) 



Pentamer 



< > 

/l-92 



FIG. 7: Equilibrium geometries of H20 n , n—3-5; interatomic 
distances in A. 



reported in Table IVTlTl 



A. Small water clusters: trimer-pentamer 

Early spectroscopic studies^ predicted the open chain 
conformation to be the most stable structure for the 
trimer. Subsequent work has suggested otherwise ; 100 ' 106 
and the cyclic trimer with C\ symmetry has been shown 
to be the more stable structure. Using the modified 
BFGS routine to perform the energy minimization, we 
found the ring structure to be slightly more stable than 
the open chain conformation with the difference in en- 
ergy being 0.861 kcal/mol (this lies within the margin of 
error for our potential fit, 0.02 eV/molcculc x 3 = 1.38 
kcal/mol.) Next, we obtained the energies and optimized 
geometries of the predicted ground-state structure of the 
tetramer and pentamer. The S4 cyclic tetramer structure 
has been shown to be most energetically favored for the 
water tetramer In this structure (Fig.[7Jc)), there are 
alternating hydrogen atoms above and below the plane of 
the tetramer ring. A puckered cyclic ring (Fig. [TJd)) has 
been predicted to be the most stable pentamer structure 
by both ab initio calculations^ and experiment^ 8 . 

Table l*VTlT1 gives the properties of the three clusters; for 
comparison, along the lines of Stern et al.^ we present 
results of select potentials along with ab initio calcula- 
tions and experiment. The notation (• • ■ ) in Tables IVIII I 
IXIII reflects the fact that all quoted distances are aver- 
aged over the cluster structure. Though our model pre- 
dicts the correct structure and energetics, the net dipole 
moment \x is once again smaller than that computed via 
other models as well as experiment, as expected based on 
our previous discussion. Scaling /i by a factor equal to the 
ratio of the experimental monomer dipole moment and 
our model's monomer dipole moment (/x„ orm = 0.6613) 
yields values that are more realistic; these arc the values 



B. Water hexamer 

For the water hexamer, it has now been established 
that there are a number of different local minima struc- 
tures that are energetically very comparable. IR spec- 
troscopic experiments on gas-phase clusters by Paul et 
a/J-iii and Liu et qZj 109 ' 110 indicate that the caged hex- 
amer structure is the most stable, while ab initio cal- 
culations have revealed that the cage, prism and book 
structures are almost degenerate, with the stability se- 
quence depending on the inclusion of zero-point energy 
differences, 117 ! 118 ! 119 ! 120 ' 121 ! 122 ' 123 Further, Tissandier et 
alr^l have used a topological enumeration technique in 
conjunction with semi-empirical PM3 methods to predict 
the global minimum energy structures. Here we examine 
the cyclic, cage, prism, chair and book structures; the re- 
sults are provided in Tables ITXl and 1X1 The prism, book 
and the cage structures are the most stable and are en- 
ergetically nearly degenerate, while the cyclic and chair 
are clearly metastable structures at K. The computed 
dipole moment \i has been scaled by [i norm . Fig. [5] shows 
the various water hexamers as obtained from our model. 
The results clearly indicate that our model is capable of 
describing the experimentally determined structures and 
relative energetics of the water hexamers. 



C. Beyond the Hexamer 

The experimental energetics and the structures of wa- 
ter clusters with six or fewer molecules have been well 
documented . 88 - 89,90 - 91,105,106,107,108,109 - 110 - 111 - 112 This is 
not true for larger water clusters (n > 10), and informa- 
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TABLE VIII: Equilibrium properties of water clusters (n=3-5). All distances in A, angles in degrees, energies in kcal/mol, 
charge in e, and dipole moment [i in debye (D). 6-31G** HF energies are given in parentheses along with select ab initio values. 





Predicted 


POL5/TZ 


POL5/QZ 


TIP4P/FQ 


TIP5P 


MCDHO 




ab initio 


Expt. 


Trimer- Cyclic 




















V 


-13.743 


-13.416 


-13.453 


-12.576 


-14.992 


-13.982 


-15.9 C 


(-17.10) 




(too) 


2.712 


2.901 


2.893 


2.912 


2.770 


2.911 




2.782 e 


2.960 h 




0.673 


1.205 


1.205 


0.417 


1.074 


1.114 




imi f 




(Zhoh) 


105.26 


















{row) 


0.949 


















(9H> 


0.221 


















Tetramer- Cyclic 




















U 


-27.308 


-25.529 


-25.665 


-23.641 


-28.431 


-27.581 


-23.8 9 


(-29.10) 




(roo) 


2.826 


2.769 


2.759 


2.809 


2.673 


2.806 




2.743 9 


2.79 h 




0.015 


0.000 


0.000 


0.000 


0.000 


0.024 




0.000^ 




(Zhoh) 


105.27 


















(roH) 


0.948 


















(«h) 


0.220 


















Pentamer— Cyclic 




















u 


-36.027 


-34.111 


-34.427 


-32.954 


-38.122 


-35.229 


-33.34 9 


(-37.70) 




(roo) 


2.850 


2.742 


2.726 


2.773 


2.657 


2.753 




2.867 9 


2.760 h 




0.634 


1.190 


1.191 


0.401 


1.219 


0.992 




0.927 ; 




(Zhoh) 


105.26 


















(roH) 


0.948 


















(9H> 


0.220 



















a Ref. [43 
''Refs [2pEl 
c Ref.[33 

Toil . 

jRefs .ll03lll04t 
gRef. llOOl. 

fe Refs. 



d Ref. 
c Ref. 



105 ld| (trimer) ;[m|UTJ (tetramer) ; [l07llT08l (pentamer). 



tion about such clusters is available mainly via classical 
potentials and quantum calculations. Hence we compare 
our results only with other computational studiesi 86 ' 96 i 97 
Maheshwary et al£^ have examined the structure and 
stability of water clusters (up to twenty-molecule clus- 
ters) using Hartree Fock as well as DFT (B3LYP) calcu- 
lations with 6-31G** and 6-31++G** basis sets; calcu- 
lations using TIP4P22 and TIP5P-2& potentials have also 
been performed for these clusters. 

The experimentally-determined^ and theoretically 
predicted^ stable heptamer conformer is a cuboid struc- 
ture with a missing corner, labeled Heptamer (a) in 
Fig. [9] This is also the lowest energy geometry as pre- 
dicted by our potential, with an unsealed dipole mo- 
ment of 1.20 D. In addition, we observe another structure 
(Heptamer (b) in Fig. [9]) to be approximately 1 kcal/mol 
higher in energy. This structure has a high dipole mo- 
ment (3.94 D), and nine hydrogen bonds, in contrast to 
the ten found in the more stable conformer. Note that the 



dipole moments reported in Table IXll for the large clus- 
ters are as obtained and have not been rescaled, since 
no experimental data is available for comparison. (We 
also would expect the deviation between theory and ex- 
periment resulting from our specific choice of atom-in- 
moleculc charge definition to "wash out" for the larger 
clusters.) The same ordering in the energies and dipole 
moments is seen in the work of Maheshwary et alf^ 

The most stable state of the water octamer in our work 
is cubic with D2d symmetry. The next most stable oc- 
tamer structure is another cubic structure with SA sym- 
metry. We observe a difference of almost 1.4 kcal/mol in 
the relative energies of the two structures; Maheshwary et 
al§£- predict the two structures to be nearly isoenergctic. 
The dipole moment is zero for both structures, with each 
structure characterized by twelve hydrogen bonds. These 
structures are shown in Fig. [9] 

The global minimum water nanomer structure can be 
described in terms of a pentamer and a tetramer ring con- 
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TABLE IX: Equilibrium properties of water hexamers I. All distances in A, angles in degrees, energy in kcal/mol, charge in e, 
and dipole moment /i in D. 6-31G** HF energies are given in parentheses along with select ab initio values. 





.Predicted 


rULo/ IZi 


FULo/CjZ 


1IF4F / tL) 


mTD rnc 

llroF 




DC 


ab initio 


Expt. 


Hexamer- Cage 




















V 


-46.497 


-41.783 


-39.297 


-45.388 


-45.388 


-43.690 


-40.76 


-45.03 / (-48.60) 




(roo) 


2.801 


2.783 


2.755 


2.863 


2.746 


2.888 




2.807 / 


2.820'' 




2.120 


2.442 


2.454 


1.788 


2.178 


2.034 




2.05 9 


1.904'' 


(Zhoh) 


105.14 


















(rem) 


0.950 


















(<2h) 


0.220 


















Hexamer- Book 




















U 


-46.492 


-42.464 


-42.771 


-40.152 


-46.680 


-43.977 


-40.38 


-44.74 ; 




(roo) 


2.788 


2.777 


2.777 


2.815 


2.688 


2.809 




2.766 / 




A 4 


2.410 


2.449 


2.430 


2.006 


2.445 










(Zhoh) 


105.22 


















(i-oh) 


0.949 


















<9H> 


0.220 


















Hexamer- Prism 




















U 


-46.465 


-41.847 


-42.135 


-39.304 


-45.805 


-44.192 


-40.97 


-45.12 / (-49.60) 




(roo) 


2.757 


2.792 


2.782 


2.819 


2.773 


2.892 




2.840 / 




/' 


2.974 


2.905 


2.931 


3.254 


2.692 


2.627 




2.701 9 





(Zhoh) 
(j"oh) 
(«h> 



104.85 
0.951 
0.219 




»Refs. ll03irLU4l . 

fe Refs. [ToglfTTol. 



nected by hydrogen bonds (Nanomer (a) in Fig. [5]). This 
structure is seen by experimental studies of Buck et al~2!l 
as well as computational studies by Maheshwary et al£& 
and Dang and Chang,— and is characterized by thirteen 
hydrogen bonds. Our potential also predicts this struc- 
ture to be the most stable. Another stationary point on 
the nanomer energy surface is the structure Nanomer (b) 
as shown in Fig. [9l This structure contains 13 hydrogen 
bonds, and can be described as a octamer cube plus a 
monomer coordinated to a corner of the cube via a hy- 
drogen bond. 

Locating the global energy minimum for larger clus- 
ters (n > 10) is a difficult task since the flat potential 
energy surface gives rise to many possible geometries 
with comparable energies. We have therefore used the 
geometries predicted by Maheshwary et al.,— TIP5P;** 
and TIP4P— (available online at the Cambridge cluster 
database website^) as starting configurations for our 
energy minimization calculations. The energies of our re- 



sulting energy-minimized structures for n = 10 — 20 (Ta- 
ble IXII[) , agree reasonably well with the calculations of 
Maheshwary et al. Rather than providing the geometries 
of all the above clusters, we have listed their important 
properties in Table IXHl the table also indicates the ini- 
tial structure that yields the minimum energy geometry 
when we perform our minimization. 

A summary comparison of our model-predicted results 
with the ab initio calculations of Maheshwary et al. is 
given in Fig. 1101 Some deviations from the ab initio re- 
sults occur at n = 3, 6, and 16. In particular, as shown 
in Table IVIII1 we underestimate the binding energy of 
the trimer, leading to the deviation in estimation of the 
incremental interaction energy at n = 3. Although the 
model is able to predict the correct ordering of the bind- 
ing energies of the various hexamers and heptamer, it is 
unable to capture the small difference in incremental in- 
teraction energy between n ~ 6 and n ~1 . However, we 
do largely reproduce the alternation in stability of the 
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TABLE X: Equilibrium properties of water cluster hexamers II. All distances in A, angles in degrees, energy in kcal/mol, charge 



in e, and dipole moment /i in D 


6-31 G** HF 

. \J rj J_\J 111 


PTlPTCnpQ arp 


envpn in nfirpnfVipQpQ alnncr with tiplppf /7ft 7757/10 valnpQ PnntnntpQ an 

tlVL.ll 111 IJtXl LlllllLoV^iJ Cl 1 W 1 1 1_. VV 1111 lj d C L Li l/ (it VVVCvJ V CXI U.CD. 1_ vJ VJ ull W Lr L> l> CLO 


in Table US 












Predicted 


rULo/ 1 Zj 


DAT C lr\ r 7 a 


llr^r/rVc^ ll.ro_r IVlUJJxlU JJL> initio rjxpt. 


Hexamer- Chair 










U 


-44.073 








(roo) 


2.846 








/' 


n ni i 

U.Ull 








(^hoh) 










(rem) 












0.220 








llcAcLIIlcl y 










U 


-43.919 


—41.875 


—42.224 


—41.368 —47.309 —44.264 —39.34 — 43.88 J 


{roo) 


2.849 


2.737 


2.720 


2.756 2.654 2.731 2.714 ; 2.756 9 




0.150 


0.017 


0.003 


0.000 0.000 0.134 0.000 


(/-hoh) 


105.25 








(rem) 


0.948 








<«h> 


0.220 









Hexamer- Cage 



<k Jr* 



ia) 



, Hexamer- Prism 

T Hexamer- Book Xjfa 



(b) 



y Jt 



Hexamer- Chair 



(d) 



(e) 



FIG. 8: Equilibrium geometries of water hexamers; inter- 
atomic distances in A. 



cluster depending on whether n is odd or even — in par- 
ticular, the enhanced stability of even n-mers relative to 
odd n-mers. 



V. DISCUSSION AND CONCLUSIONS 

We have presented a new dynamical CT-EAM poten- 
tial for modeling water and its polymorphs at the atomic 
level. We have based our parameterization on ab ini- 
tio data; in particular, atomic charge fluctuations have 
been modelled with reference to the local chemical en- 
vironment, using Lowdin population analysis to repre- 
sent atomic charges. Depending on its immediate coor- 
dination environment, each atom is assigned to a cluster, 




^ Oc la m or- 



(d) 




•v. 

t *, /l .96 

(c) 

V 



r 



FIG. 9: Equilibrium geometries of H^On, n=7-9; interatomic 
distances in A. 



with this identification being crucial to our formulations. 
Cluster identification is effected via a radial cutoff, and 
total cluster charge is based on the size of the cluster and 
relative positions of neighboring clusters. This charge is 
in turn partitioned among the constituent atoms. 

Our technique is sufficiently flexible to account for very 
different charge states of clusters and individual atoms. 
The radial cutoff chosen to define our clusters, r cut = 1.5 
A, is significantly larger than the O— H equilibrium dis- 
tance in the monomer and dimer (~ 0.95 A.) Conse- 
quently, the model is easily capable of describing non- 
perturbative charge transfer. 

We note that a number of important effects have been 
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TABLE XI: Equilibrium properties of water clusters for n~7- 
9. All distances in A, angles in degrees, energy in kcal/mol, 
charge in e, and dipole moment fj, in D. 





Predicted 


TIP4P" 


TIP5P 6 


ab initio c 


Heptamer (a) 










u 


-58.259 


-58.271 


-57.910 


-60.53 


(fnn) 


2.802 


2.762 


2.738 


2.884 




1.20 






1.35 


{ ^- hoh) 


105.09 






106.26 




0.951 






0.950 




0.219 








Octamer Z>2d 












-74.325 


-73.090 


-72.535 


-76.01 


(rnn) 


2.822 


2.746 


2.712 


2.877 




0.00 






0.00 


(^hoh) 


105.13 






106.482 


(von) 


0.951 






0.951 


(<7h) 


0.219 








Nanomer (a) 












-84.124 


-82.401 


-83.622 


-85.05 


{roo) 


2.842 


2.741 


2.696 


2.869 




0.99 






1.69 


(Zhoh) 


105.16 






106.39 




0.950 






0.950 


<9H> 


0.219 









a Ref.[§3 

<>Ref. 
c Ref. 




10 12 

Cluster size n 

FIG. 10: Incremental interaction energies of water clusters, 
AE = E„+i — E n — Ei, as a function of cluster size n. Note 
that AE = Un+i — U„, where U n is the binding energy for a 
cluster of size n. UHF results from Ref. 18611. 



omitted in this initial implementation. This was done 
in order to focus attention on the physics of the charge- 
transfer EAM model framework itself, rather than the re- 
finement of a model water potential per se. For example, 
the current parameterization does not yet impose the cor- 
rect asymptotic dissociation behavior on cluster subsys- 
tems, which we have argued is critical to a proper descrip- 
tion of reactive dynamicsi 69 i 71 A related issue concerns 
the omission of several ionic species, believed to be impor- 
tant in defining the hydrogen network in water, from our 
paramcterizations: these include the Zundcl (H5OJ) and 
Eigen (HgO^ - ) cations^ A final important simplification 
concerns the use of atom-in-molecule charges as prox- 
ies for the shape function modeling of the AIM charge- 
density distributions. It is clear that further work taking 
account of these various factors will be necessary in or- 
der to successfully study complex kinetic processes such 
as those involved in ion solvation, enzyme catalysis, and 
proton transport. This work is presently underway. Ad- 
ditionally, it should be noted that our approach does not 
incorporate a quantum mechanical treatment of the ac- 
tual electron or proton transfer processes.*^ 

Notwithstanding the simplicity of this initial model, in 
tests on small water clusters, our results agree very well 
with experimental and ab initio data. Importantly, our 
model captures the transition from planar ring pentamer 
structures to three-dimensional complex hexamer struc- 
tures, an essential structural test for any successful water 
potential. In this context, it is worth noting that an envi- 
ronment dependent dynamic charge potential motivated 
by the present work has also been developed recently 
for silica. This potential successfully matches ab initio 
results in its ability to predict the ground-state energy, 
geometry and failure mechanisms of silica clusters.*^ 

More generally, this work represents a successful ap- 
plication of many-body embedded atom concepts to the 
modeling of a highly polarizablc molecular system, and 
thus a significant departure from traditional approaches 
to developing water potentials. It is remarkable that even 
this relatively simple implementation of CT-EAM repro- 
duces cluster structures and energetics consistent with 
the best previous potentials, while providing a theoret- 
ical roadmap for implementing true charge-transfer dy- 
namics. This ability of an embedded-atom approach — 
originally designed for describing many-body effects in 
bulk fee metals — to model the structure of a molecular 
system can be understood as a direct consequence the 
CT-EAM framework's underlying density functional con- 
struction. DFT, with its emphasis on electron densities 
as the fundamental variables of the theory, acts as a mul- 
tiscale mechanism for incorporating quantum mechanical 
bonding effects and excitations within a nominally clas- 
sical potential. 

We believe that the unique combination of features de- 
scribed here will ultimately enable CT-EAM potentials 
to successfully capture many-body and electrostatic ef- 
fects, in both static and dynamic contexts, for a wide 
variety of biophysical and materials systems, including 
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TABLE XII: Predicted equilibrium cluster-averaged properties of water clusters for n > 10. All distances in A, angles in 
degrees, energy in kcal/mol, charge in e, and dipole moment fi in D. "Geometry" refers to the starting geometry for the energy 
minimization, as described in the text. 



n 


u 


(roo) 


A« 


(Zhoh) 


(ran) 


(m) 


Geometry 


10 


-96.311 


2.841 


1.81 


105.12 


0.951 


0.219 


TIP5P 


11 


-105.098 


2.849 


2.54 


104.97 


0.952 


0.219 


TIP5P 


12 


-119.291 


2.828 


0.00 


104.87 


0.953 


0.219 


Ref. [86] 


13 


-126.951 


2.835 


1.65 


105.11 


0.951 


0.219 


TIP4P 


11 


-144.868 


2.835 


1.86 


105.01 


0.952 


0.218 


TIP4P 


15 


-152.589 


2.850 


1.98 


105.12 


0.951 


0.218 


TIP5P 


16 


-162.153 


2.829 


0.00 


104.73 


0.954 


0.219 


TIP4P 


17 


-174.362 


2.848 


3.03 


105.02 


0.952 


0.219 


TIP5P 


18 


-192.442 


2.836 


1.85 


104.93 


0.953 


0.218 


TIP4P 


19 


-201.608 


2.842 


2.97 


105.03 


0.952 


0.218 


TIP5P 


20 


-214.889 


2.833 


0.17 


104.93 


0.953 


0.218 


TIP4P 



nanoscalc systems possessing mixed molecular and bulk 
features. Future applications will include studies of the 
crystalline polymorphs of water and the thermodynamic 
and structural properties of the liquid, as well as investi- 
gations of dynamical processes such as ion solvation and 
proton transport. 
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